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ABSTIACT - la LABDSAT laage ry, spectral and 
spatial information oaa be uaad to datoat 
tba drainage aatvorlc aa wall aa tha rela- 
tive elevation nodal la aottatalaoaa tar- 
rain. To do this, alzad information of no- 
tarial raflaotaaoa aad topocraphio 
modulatioa la tba original LAIOSAT imagery 
east ba aaparatad first. Froa tha aatarlal 
raflactaaea Information, big visible rlvara 
caa ba dataotad. Froa tha topocraphio mo- 
dulation information, rldcas aad valleys 
can ba dataotad aad aaalcaad relative ele- 
vations. A ooaplata elevation nodal oaa ba 
generated by latarpolatiac values for aoa- 
rldga and non- valley plzals. Tha saall 
straaas aot dataotabla froa aatarlal re- 
flaetaooa information caa ba loeatad la 
tha vallaya with flow dlractlon known tie x 
tha elevation nodal. Finally, tha flow di- 
rections of blc visible rivers oaa ba In- 
ferred by aolviac a coaslstaat labellac 
problea based oa a set of spatial reasonlnc 
coastralats. 


X. latradacUoa 

. It is a coaaon taste for a photointer- 
preter to ezaaine the spatial pattern oa aa 
aerial laace and by appropriate interpreta- 
tion be able to tall the elevation of oaa 
area relative to another aad ba able to in- 
fer the streaa network aad the dralaace 
network even thou(h soaa of tha straaas aay 
be below tha resolution of tha sensor. 
There is a wealth of laforaatloa la spatial 
patterns oa aerial iaagery but most compu- 
ter data proeessiac of reaotely sensed Im- 
agery, being Halted to plzel spaotral 
characteristics, does not aake use of It. 

In this paper, we describe a prooadure 
by which the streaa network aad relative 
elevation model can be laferel froa a 
LAMDSAT scene of mountainous and hilly ter- 
rain. The processing has a nuaber of dis- 
tinctly different steps. First to appro- 
priately prepare the iaagery for processing 
ve must destripe it aad perfvra haze remo- 
val. Cestrlplng caa be done by the Kora 
and Voodhaa [1979J technique. Base removal 
can be done by the Switzer, Kowallk aad 
Lyon technique. These two steps 
constitute the preprocessing and are not 


dlsoassod la this paper. 

To a first order effect, after prepro- 
cessing the cause of the intensity value at 
aay plzel Is due to the combined effect of 
the aaclo at which the sua illuminates the 
ground patch corresponding to the plzel aad 
the reflectance of the surface material oa 
the ground patch. To aake seaae or the 
spatial pattern first requires separatlag 
these two effoots. For this purpose the 
Ellaaoa, Soderbloa aad Chaves [1991] tech- 
nique caa be used to create two iaages froa 
the oae LAIOSAT laage. The f irst laage . Is 
a reflectance laage aad the second laage'' is 
a topographic aodulatloa laage and has la- 
foraatloa related to surface slope aad sua 
illumination. The details are given la 
Section 2. 

As discussed la Section 3, tha reflec- 
tance laage caa be used by the Alfoldl aad 
Munday Cl 978] procedure for identification 
of all areas of water. The topographic mo- 
dulation laage oaa be used to Identify the 
ridges aad the valleys. This is discussed 
la Section A. kith the valleys identified, 
each valley plzel aay be assigned a rela- 
tive elevation which laoreases aa the val- 
ley path froa the plzel to the river it 
eaptles la increases. Ildges ' aust be as- 
signed elevatloaa higher than their neigh- 
boring valleys aad each ridge plzel caa be 
assigned a relative elevation which la- 
creases oa the ridge path froa the plzel .to 
the saddle polat where the ridge orosaaa a 
valley increases. The ridge valley eleva- 
tion assignment procedure la discussed la 
Section S. Once ridges aad valleys have 
been located and assigned relative eleva- 
tions, a complete elevation model can be 
generated by interpolating values for aoa- 
rldge and non-valley plzels. The Interpo- 
lation procedures are discussed la Seotlon 
6. The final eleaent of the spatial rea- 
soning is the asslgnaent of streaa or river 
flow direction for those water bodies which 
were directly identified by reflectance 
properties. This Is discussed in Seotlon 7 
and 8. In the remainder of this section, 
we review previous work done by saptial 
reasoning Investigators. 
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Sine* the launch of tha first Earth 
Resources and Technology Satalllts (CRTS, 
later renaaed LANDS AT) la July 1972, auoh 
work la reaote sensing has been done by us- 
ing pattern analysis and picture processing 
techniques for laage classification, re- 
storation and enhaaceaent. Pew people have 
►rled the scene analysis or artificial In- 
telligence approach to deserlbe the laage 
in teras of the properties of objects or 
regions la the laage and the relationships 
between then. Ehrlch Cl 977 1 found global 
llneaaents by partitioning the laage into 
windows and applying long, straight linear 
filters at different orientations in each 
window to extract local evidence. Dynaalo 
progressing (Montsuarl, 1971; nartelli, 
1972] was then used to fora ooaplete global 
llneaaents. TanderBrug (1976] tested vari- 
ous detectors to get linear features la sa- 
tellite laager?. This was only at the lo- 
cal level. Later TanderBrug (1977a] used 
relaxation to reduce noise la the output. 
Finally TanderBrug (1977b] defined a aerlt 
function that can be used to select pairs 
of segaents to be aerged so that local line 
detector responses can be linked together 
into a global representation of the curves. 
His work Is closely related to the Shlral 
(1973] technique which eaployed sequential 
line following to find edges in soeaes con- 
taining polyhedra. LI and Fu (1976] used 
tree graaaars to locate highways and rivers 
froa LANDSAT pictures. The above Investi- 
gations deal with the extraction of all the 
linear features froa an laage, but they do 
not deal with the Interpretation of these 
linear features. In the following Investi- 
gations, knowledge about the desired fea- 
tures are considered crucial la such ana- 
lyses. 

Bajcsy and Tavakoll (1975] argued that 
an laage filter is not aeaalngful unless 
one has a world aodel, a description of the 
world one Is dealing with. They recognized 
objects notching this description and fil- 
tered then out. This strategy Is used to 
sequence the recognition of bridges, riv- 
ers, lakes, and Is) ds froa satellite pic- 
tures. Nagao and Ma^ luyaaa (I960] built an 
inage understanding systea that autoaati- 
cally located a variety of objects in an 
aerial photograph by using diverse know- 
ledge of the world. It Is one of the first 
laage understanding systeas that has Incor- 
porated very sophisoated artificial intel- 
ligence techniques Into the analysis of 
coaplex aerial photographs. Plsohler, Ten- 
enbaua and Wolf [1961] designed a low-reso- 
lution road tracking (LIRT) algorltha for 
aerial laagery. The approach was based on 
a new paradlga for coablninr local iaforaa- 
tlon froa aultiple sources, asp knowledge, 
and generic knowledge about roads. The fi- 
nal interpretation of the scene was ac- 
hieved by using either graph search or dy- 
namic progressing. 


Slallarly, knowledge Is lsportant In 
our problea whloh requires analysis both at 
the local and global levels. Looal level 
analysis will be discussed la Seotlon 2 to 
A; global level analysis will be dlsoussed 
in Seotlon 5 to 8. 


i- IllUNlaaUBC nodal 

The brightness and darkness la eaoh 
band of LANDSAT laage s coae froa two aaln 
sources. First, they oan be due to aaterlal 
properties. For exasple, in the spectral 
region (.6 - 1.1 pa) of band 7, water bod- 
ies absorb infrared radiation, so they ap- 
pear as clearly delineated dark bodies; 
living vegetation reflects strongly in this 
portion of the Infrared, so areas of living 
green vegetation appear as bright regions. 
Second, they any be due to topography and 
sun lllualnatlon angle effects. The aoun- 
taln side facing to the sun appears es a 
bright region; the aountaie side facing 
away froa the sun any appear as a shadow or 
dark region. Unfortunately, the LANDSAT 
data values are soae coablnation of these 
two effects. Eliason, Soderbloa, and Cha- 
ves (1981] address this problea by defining 
an lllualnatlon aodel. In the following, 
their general theory about the brightness 
in LANDSAT laagery will be Introduced, and 
extraction of aaterlal reflectance and to- 
pographic aodulatlon inforaatlon based cn 
clustering on ratio iaages will be de- 
scribed. 

The original LANDSAT laage 8* aeasur- 
iog the aaount of light reflected froa a 
surface at pixel (r, c) for wavelength w , 
1 a A, 5. 6, 7 Is : 

8* ( r, e, w ) • I(r, e, w.) • Tp(r, c, 

p) *. H (w ) 1 1 

where R is the brightness of the scene if 
Che surface were flat, Tp Is the aodulatlon 
of the brightness introduced by topography, 
and H Is the haze due to ataospherlc scat- 
tering. The laage of R Is called aaterl&l 
reflectance laage, and the inage of Tp li’ 
called topographic aodulatlon laage. p Is 
the photoaetrlc function which depends on 
the phase angle, the Incidence angle, and 
the angle of eaergence (vildey, 1975], but 
which does not depend on w.. Two assusp- 
tlons are aade here: 

1. Tp Is independent of aaterlal proper- 
ties and wavelength. 

2. The photoaetrlc function p is indepen- 
dent of wavelength. 

After H is calculated by the Switzer, 
Eowallk and Lyon [1981] technique, for each 
band, H(w ) is subtracted froa 8'(r, c, w.) 
at ail pixels to get the haze-corrected la- 
age 

B(r, c, w.) « B’(r, e, w.) - H(w.) 

a R(r, c, w t )*» Tp( r , c, p). 
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Because Tp Is Independent of 

wavelength. In the ratio lnage of two bands 
with wavelength w^ and Wj , each pixel (r, 
c) haa gray level 

B(r, c, Wj) / B (r, c, w 2 ) 

« (R(r, c, w ) • Tp(r, o, p ) ) / (R(r, c, 
w 2 ) • Tp ( r , o, p)) 

■ R(r, c, w - ) / R (r, c, w_) , 
which Is Independent of Tp. Thus, the to- 
pographic information Is reaoved. 1 simple 
demonstration of this theory Is that, In 
the ratio Image after all diffuse lighting 
has been removed, all shadows disappear. 

Ellasoa, Soderblom, and Chavez estimate 
the material reflectance Image R for each 
w by clustering using different ratio Im- 
ages as features. The result In a cluster 
Image Cl (r, c). For each w the average 
brightness value of all the pixels In each 
cluster is taken to represent R for their 
respective cluster. Two basic assumptions 
are Inherent: 

1. All image elements that group in a ra- 
tio cluster represent a single material. 

2. The topographic slopes of all elements 
in a cluster are symmetrically distributed 
toward and away from the sun, such Irat 
their average brightness can be used to es- 
timate the brightness of that mat*rl«l on a 
flat surface. 

The material reflectance Images will be 
used In detecting visible rivers described 
In Section 3. 

From the cluster image, four material 
reflectance images R(r, c, w ) can be esti- 
mated because, for each cluster, four aver- 
age brightness values can be calculated 
from the four bands, l.e., for 1 * 4, 5, 6, 
7, if Cl (r, c) » k at pixel (r, c), then 

EtZy B ( x , y , w ) 

Cl( x, y ) >k l 

R(r,c,w ) * — — 

Cl(x,y)«k 

For any pair of the haze-corrected image 
3(r, c, w^) and material reflectance image 
R(r, c, w.) , 1 * 8, 5, 6, 7, the topograph- 
ic modulation image Tp can be calculated 
simply by taking the ratio of B over R. 
This topographic modulation image is an im- 
age whose tonal variation is unambiguously 
identified with surface slope and sun illu- 
mination angle. 


1. Patecllflfl al Visible Rivers 

Because one pixel in the LANDS AT image 
represents approximately a 57 meter by 80 
meter area on the ground, the resolution is 
low. For the most part, it is not possible 
to directly observe the drainage network 


of the LANDSAT data. If rivers or lakes 
are visible to the humans, they can be de- 
tected by spectral Information as described 
In this section. Ve call such rivers (in- 
cluding big lakes) ‘visible rivers* On the 
other hand. If the streams are not risible, 
they can only be defeated by spatial Infor- 
mation and we call them ‘invisible streets/ 
For any window over the LANDSAT Image, If 
one can deteat some visible rivers by using 
speotral information and detect ridges and 
valleys by spatial pattern, then It Is pos- 
sible to continue to look for invisible 
streams by using spatial information. Ex- 
amples of visible rivers and Invisible 
screams are shown In Figure 1 . The image 

was taken in April, 1976 over areas in Ni- 
cholas County, U . Fa. and neighboring 
counties . 



Figure 1 - LANDSAT scene in W. Va. 
a indicate a visible river, and 
b, c indicate invisible streams 


Once the material reflectance image is 
created by the technique In last section, 
It can be used to identify visible rivers. 
In the spectral regloa (.8 - 1.1 pm) of 
band 7, water bodies absorb infrared radia- 
tion, so visible rivers appear as dark 
curves, and lakes appear as dark regions. 
In the material reflectance image of band 
7, these dark features become more clear 
because shadows are removed. However, not 
all dark features are water bodies; the 
real water oodles can be identified by the 
following process [Alfoldi and Kunday, 
1 978]. 

(1) A band 4 green coefficient x of 
every pixel is calculated as the ratio of 
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the radiance of band 4 over the radiance 
sun of bands 4, 5 and 6. Similarly a band 
5 red coefficient y is calculated for every 
pixel. x and y are called LANDSAT cbromat- 
icity coordinates. 

(2) In this aoordinate system, Hunday 
[1974] has determined a curve (Figure 2) 
which is the loous of the positions of 
chromatlclty values of water bodies. If, 
for some pixels, the z, y values calculated 
in 1 are close to this curve, then those 
pixels can be identified aa portions of wa- 
ter bodies. 



Figure 2 - Chromatlclty plot 


JL. Rldge-7alley &ju AXLSl Invisible Stfsaaa 

After visible rivers are detected by 
both spectral and 20 spatial Information, 
invisible streams can be detected by 30 
spatial information. In this nountainous 
area, water flows through valleys, so that 
the drainage network of invisible streams 
is a substructure of the valley network, 
therefore, the approach suggested here is 
to first get a relative elevation model, 
then extract the valley network from this 


elevation model, and finally extract the 
drainage network of invisible streams from 
the valley network. In this section, we 
desorlbe how to extract shadowed and bright 
areas, create linear features on the bord- 
ers between these areas, and classify these 
linear features into ridge and valley seg- 
ments. In the next two sections, we dis- 
cuss how to generate a relative elevation 
model. The extraotlon of valley network 
and the network of invisible streams will 
be discussed in Section 9. 

In the topographic modulation image, 
bright areas Indicate that the surfaces are 
facing to the sun, dark areas Indicate sha- 
dows. In order to detect valleys and ridg- 
es, it is necessary to first segment the 
image into regions of shadowed and bright 
areas because valleys and ridges exist on 
the borders between these regions. Gray 
level thresholding can be used to determine 
shadowed and bright areas. We use the Wa- 
tanabe [1974] technique in a recursive way 
to select thresholds. The details are giv- 
en in Wang and Hers! 1c k [19*2]. The con- 
nected components of them are shown in Fig- 
ure 3 • 



Figure 3 - a. Connected components of dark 
regions 
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Figure 3 - b. Connected components of 

bright regions 


Kext, tho perimeters of these bright 
and shadowed regions are segmented Into 
"border segments" according to their left re- 
gions, right regions, and orientations. A 
border segment is a maximally long sequence 
of connected pixels which are on the border 
between two given regions. Because the de- 
tection of ridges and valleys Is highly or- 
ientation-deoendent and the sun illumina- 
tion comes from east in Figure 1 , each 
border segment is further broken into sev- 
eral pieces according to orientation: all 
the east-west parts can be separated from 
the north-south parts. The result Is shown 
in Figure A . 



Figure A -Border Segments 


As the sun illumination comes from 
east, those border segments which are val- 
ley segmeats or ridge segments can be iden- 
tified according to the brightness of their 
left and right regions. 8ecause most of 
the trees in this area in April are unfoli- 
a”ed, tha strongest region boundaries are 
shadow bou darles rather than tonal boun- 
daries, and the strongest boundaries are 
those at the extremes of steep slopes 
oriented normal to the sun direction. Be- 
cause the sun illumination is predominantly 
east-west, a boundary that is dark on the 
left and bright on the right will corres- 
pond to a ridge, and the reverse will cor- 
respond to a valley. 

For east-west region boundaries, tha 
above ridge-valley inference mechanism 
fails. Where east-west boundaries exist, 
some are ridges and some are valleys. To 
classify these east-west border segments 
correctly, it requires elevation informa- 
tion. As shown In Figure 5, if end a of 
the valley segment VI Is higher than end b, 
Z can be determined to be a valley. Also, 
if end a is lower than end b or about the 
same, Z can be determined to be a ridge. 
The results of ridge-valley finding are 
shown in Figure 6. Assignment of relative 
elevation to ridge and valley is discussed 
in the next section. 


C 
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Figure S - Classifying aaat-wast bordar 
ii|itBta. ft ia a vallay segaeat; K1 la a 
rldga segaent 




bordar aagaaata which ara identified aa 
ridgaa. 


a- Btlitlrt IliritiflM aL Jltfu aa 4 f hilar 
Sagaaata 

Tha dataetloa of tha rldga aad valley 
aagaaata aa dlaeuaaad la tha laat aaotioa 
oaljr aaalgaa a rldga or valley labal to 
thaa aad doaa aot aaalga relative eleva- 
tloaa to thaa. Ia thla aaotioa, wa de- 
aerlba how to aatlaata thalr ralatlva ele- 
vatloaa. Plrat wa will daaoriba a aodal 
which can do tha alavatloa aaalgaaaat Job, 
than wa will glva tha equations of eleva- 
tlon aaalgaaaat. 


Figure 6 - a. Vallay aap oonaistlag of tha 
bordar sagaants which ara idaatlflad aa 
vallaya. 


laauaiag that wa hava a atraaa network 
la a aouatainoua araa, aad wa .know whara 
tha blggaat rlvara ara, wa oaa traca tha 
aatwork, atartlag froa tha blggaat rlvara, 
to find tha flow dlraotloaa of all tha 
atraaa aagaaata baeauaa watar always flows 
froa higher locations to lower locations.. 
Ia other words. If the vallay sagaaota de- 
tected la the laat aaotioa foraad a net- 
work, than starting froa tha visible rlvara 
detected In Section 3, wa eaa traca the 
network aad assign relative elevations to 
all tha aagaaata. Unfortunately, tha ob- 
served vallay sagaants do not fora a net- 
work; there ara aany gaps. la shown la 
Flgura 7, if It ia dark on tha right and 
bright on tha left of atraaa Yb, than Vg 
oaanot be detected due to tha shadow on tha 
right of Vb, and a gap exists between Yb 
aad a saaller atraaa Vs. 




Vb 



Figure 7 - The gup between a saaller and a 
larger aereaa 


The knowledge that the cross-sections 
of Teller* are V-shaped eaa be need to 
bridge tho gape.- If one look* at topo- 
graphic aapa, the elevation oontoure of 
valleys each as ia Figure 8 can bo fre- 
quently found. Thus, if one drawn a line 
ab perpendicular to the Talley Fa, the de- 
rations are increasing froa point o to 
point a , and also froa point o to point b. 
However, if a ridge point is encountered 
during the process, the iaor*«slng has to 
atop because the elevation starts to de- 
crease. Thus the route of growth is di- 
rected both by the valleys and by the ridg- 
es, in other words, by global iaforaation. 



Figure 8 - The elevation pattern of valleys 
and its relation to elevation growing 


Applying this idea to Figure 7 and assuaing 
that growing propagates away froa valley 
sagaent Vb, the end a of valley segaent Vs 
will be touched first by this growing, and 
it is deduoed that end b of Vs aunt be 
higher than end a. This ia the basic idea 
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for deteralning the highar-lowcr ends of 
all the valley segaeats. The elevations of 
all the points la oao segaent «aa bo calcu- 
lated if we know its elope. On the other 
hand, rldgeo get elevatleaa when the crow- 
ing stops at thee. How, wo will give the 
staple aquations of elevation asslgnaeat. 


Our elevation growing aodel slaply as- 
sume that elevatloa laoreases aoaotonlcal- 
ly froa valleys to ridges or along valloy 
segaeats froa rivers to the saddles where a 
valley orosses a ridge. It can be used for 
asslgalng initial relative elevatloaa to 
each pixel. Because no atteapt is aade to 
realistically account for the topographic 
shape of th. hillsides froa the valley to 
the ridge, the initial relative elevations 
will be aore accurate for the ridge or val- 
ley labeled pixels than the non-ridge and 
non-valley labeled pixels. Section 6 dis- 
cusses a aore realistic procedure for hill- 
side elevatloa estlaatlon using the ridge 
valley elevations calculated in this sec- 
tion. 

There are. two ways a pixel eaa get as- 
signed an elevation depoadlng on whether 
the pixel belongs to a valley segaent or 
whether the pixel does not belong to a val- 
ley segaent. Let IT be the set of valley 
segaeats. Two slopes are assoolated with 
each valley segaent Vs ia 0: Sv(Vs) and 

Sp(Vs). Sv(Vs) is the slope along Vs it- 
self. Sp(Vs) is th* slop* of lines outside 
of Vs aid perpendicular to Vs. 

Th* elevation growing aodel constructs 
the elevation function El; Zr Z Zc -> Zp, 
where Zr is the set of row coordinates, Zo 
is the set of ooluan coordinates, and Ipais 
th* set of soro and positive integers. If p 
is a pixel belonging to a valley segaent Vs 
and pi is .the lower end pixel Identified as 
in Figure 7, then 

El(p) • El(pl) ♦ 3v(Vs) • OlstCp, pi) 
where Dint is the Euclidean distance bet- 
ween two pixels. 

If p does not belong to any valley seg- 
aest, and its elevation is originated froa 
pixel pr of valley segaent Vs, then 

Fl(p) • EH pr ) « Sp(Vs) • Dlst(p, pr). 

In a snail area, on* can assua* th* 

elevations of visible rivers are lowest. 
Assigning aoa* initial elevatloa values to 
the pixels of th* valley segaenta classi- 
fied as visible rivers, the elevations of 
all the other pixels in the laage window 
can be related to th* initial elevations of 
visible river segaenta by repeatedly using 
the above two equations. The relative 
heights of valley segaeats created by ele- 
vation growing aodel are indicated by ar- 
rows in Figure 9, and th* ground truth is 
shown in Figure to. 
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Figure 9 - Relative elevations of valley 
segments. The arrow are fro* high ends to 
low ends. 



Figure 10 - Stress asp (treated froa ground 
truth. 


*elf la formed by the junction of several 
rltfgea that radlat* outward from the peak. 
(The idealised situation represented la 
Figure 11 ahewa four symetrltally oriented 
ridges; la our area, real peaks are often 
formed by Juaotions of two or three ridg- 
es.) Rldgea of oourse are separated by 
valleys, so the higher tips of valley seg- 
ments tend to point toward peaks. The 
ridge segments laterseet to fora a peak, 
whereas valley segments tend to point to- 
wards peaks, without aetually joining. la 
thin subsection, we dlsouaa the criteria 
which can be used to identify peak junc- 
tions. 

Beoause ridge segments are the major 
features of peaks, we make the constraint 
that the number of ridge segments at a 
junction is larger then the number of val- 
ley segments. For many situations, It 
seems reasonable to relate the heights of 
peaks to the lengths of ridges that form 
the peaks. For our class of topographic 
forms (for example), It is unlikely that 
very high peaks can be fermed by the Inter- 
section of very short ridges. As a result, 
to exclude very low peaks and false peaks 
from consideration, we impose a rather ar- 
bitrary constraint opou definitions of 
peaks. Currently, we define a peak Junc- 
tion as a Juaotlon composed of four border 
segments, with the number of its ridge seg- 
ments larger than the number of valley seg- 
ments, and the length of its longest ridge 
segment longer than 800 meters. The peaks 
thus located in Figure f are narked as tri- 
angles In Figure 6.b. The eorrespondenoe 
between this result and the topographical 
map is suprlalagly good. 



1.1 td&atlflsiuon aX. £uk JuacUaai 

Figure 11 - Idealized relationships between 
When several valleys and ridges point peaks, valleys, ridges, 
toward a Junction, very often this Juatlou 
is a peak (peak at Junction). The peak it- 
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la tho laat section all pixels vara 11 - 
sigaed alavafeloaa, bat because realistic 
shape of the hillsides fro* valleys to 
ridgsa vara aot taken late account, only 
eh* ralatlr* elevation* of the ridges mod 
valleys ara hold to be aoourata. Za tbla 
aaotloa we daaorlba a few Interpolation 
procedures which peralt aora realistic ele- 
ratloa asslgaaeut to non-valley aad aon- 
ridge plxala. 

Tha flrat interpolating aurfaea baa tba 
given elevation values at ridges aad val- 
leys aad baa a 3 X 3 digital Laplaolaa of 
zopo at all noa-rldg* aad non-valley plx- 
ala. Tbla trill b* referred to aa tba La- 
plaelaa aurfaea. Tba systea of llaear 
equations which tbla ooaatralafe give* rlao 
to eaa b* written aa 

A x * b. 

Tba tractor x la tha aolutloo and represents 
the valuaa to ba aaalgaad to aaeb * varia- 
ble * (aoa-rldg* non-valley) pixel la tba 
alavatloo nodal. Tba A aacrlx la daflaad 
by applying tba digital Laplaolaa aaatc op- 
erator (Figure 12) to aaeb variable pixel. 
A aaalc operator la applied to a pixel by 
placing tha aaalc over tba laaga ao that the 
central (large poaltlve) aaak value la di- 
rectly over tha pixel whoa* value la to ba 
conputad. Tba pixel value la changed to 
■aka the «ua of tba aaak value* tlaaa tba 
corresponding laaga valuaa uadar the* equal 
to zero. For tba Laplaolaa surface only, 
Meuaann boundary conditions ar* enforced 
along tba outside rows aad eoluans of the 
elsvatloa aodel iaage. That la. the outer- 
most row or coluaa is repeated ao that the 
■ask operator can be applied to the outside 
pixels. There is one row in A for saeh va- 
riable pixel in the elevation aodel aad oae 
coefficient valur la that row for each va- 
riable. A la a sparse aatrlx alao* ao va- 
riable la constrained by aore than four 
other variables (due to the definition of 
the digital Laplaolaa aaak operator). The 
b vector la tha right band side of each of 
the linear equations la the aystea. The 
constants on the left head side of each 
equation (that result fro* applying the La- 
placlan operator to a variable pixel tbit 
has a kaowa pixel t-adjaoeat to It) are 
carried to the right bead aid* ard appear 
la b. For equation* representing variable 
pixels not t-sdjaoeat to known pixels, the 
corresponding b eleaent is zero. 


• 1 

-t « -1 
-t 

Figure 12 - A digital Laplaolaa aaak 


The second interpolating aurfaea baa 
tbs glvaa boundary values aad alalalzaa tba 
quadratic variation of tba resulting sur- 
face iCrlaaoa, 1961]. The boundary condi- 
tions with wbl04 the surface aust agree are 
depth valuaa along the zero-crossings. l'f 
the surface elevation function la C aad 
aubaorlpts denote partial differentiation, 
than the final surface E nlnialzes 

♦ «%,) «. V 

Sinoe the surface function oaa be converted 
to a discrete grid forest, the differential 
operators eaa be converted to difference 
operators, aad tba double integral caa ba 
oonvarted to double auaaatioa, the solution 
of the above fuactioa eaa be foraed by set- 
ting up a dlaorete corresponding set of li- 
near equations 

Q x e b. 

The x and b vectors have the saae aeanlng 
as la the Laplaolaa case and are construct- 
ed slallarly. The Q aatrlx Is likewise si- 
milar to the A aatrlx of the Laplaolaa. 
Instead of using Veuaaaa boundary condi- 
tions at the edge of the iaage, the quad- 
ratic variation surface la defined by using 
special peaks to fit the rows and eoluans 
near the outside edges. The six aasks 
(Figure 13) are rotated as necessary and 
applied to the only appropriate variable 
pixels .of the elevation iaage to define Q. 
Mask two la applied to corner pixels, aaak 
three Is applied to pixels la the outside 
row or eoluan that are adjaoent to a corner 
pixel, aaak four is applied to other pixels 
In the outside rows aad eoluans, aaak five 
Is applied to pixels la the next-to-the 
outside row and eoluans tbat are 6-adJaceat 
to corner pixels, aask six Is applied to 
other pixels in the next to the outside 
rows and eoluans, ‘ad aask 1 is applied to 
all other variable pixels In the Iaage. 
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Figure 14a. Elevation Modal by Method 1 
Laplaclan Maalc 


Figure 13 - Six saaks for tbe quadratic 
variation aethod. 


The third kind of interpolation surfao- 
ea can be created without using any aaak. 
For each non-boundary pixel, we can flrat 
find ita distances to the nearest valley 
pixels and nearest ridge pixels. Froa 
these distances and the elevations at these 
nearest valley pixel and nearest ridge pix- 
el, either a linear, cubic, or fifth c der 
fit interpolation can be used to ca'^ulate 
the elevation of this non-boundary pixel. 
If cubic fit is used, the first order deri- 
vative is zero at ridge and valley pixels. 
If fifth order fit is used, both the first 
and second order derivatives are zero at 
ridge and valley pixels. The resulting 
surface plots of these elevations are shown 
in Figure 14. 



Figure 14b. Elevation Model by Method 
Quadratic variation 
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Flguro IMe. 
Linenr fit 


Elevation Model by Method j, 


Figure Ite. Elevation Model by Method 3i 
Fifth order fit 


Thus far, we have found the relative 
elevations of valley segments and generated 
several Interpolation surfaces. More about 
Invisible streams Is discussed In Seotlon 
9. In the next two section, we discuss the 
assignment of flow directions to visible 
rivers detected In Section 3. 


£. Flow Direct lpna of Visible Rivers 
and Co nstraints a t Junctions 

a* 

In Section 5, an elevation growing mo- 
del was used to find relative elevations of 
Invisible streams In valleys. It remains 
to find the now directions of visible riv- 
ers which are assigned constant elevations 
In the elevation growing model. This prob- 
lem of assigning labels of (upstream, down- 
stream) to the visible river segments 1.' 
very much like the Waltz [1975] problem of 
labeling edges of polyhedra objects, and we 
need to find constraints applicable to 
streams . 


Figure 1 Ud . 
Cubic fit 


Clevatlon Model by Method 3. 


It Is believed that several stream seg- 
ments Joint at a Junction with certain or- 
ientation and length patterns. The most 
obvious and Important one Is the configura- 
tion similar to Figure 7. It is plotted as 
Figure 15. When a smaller stream S2 flows 
Into a larger stream SI S3, vary often the 
angle between S2 and SI is less than 90 de- 
grees. General rules about flow directions 
at Junctions are given in Table 1 . 
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Two ooastpslnt relations osa bo forual- 
If stattid on the basis of Tabl* t as fol- 
lows. Ono is about all 2-6upl*a 0 f sfepoaa 
seguents chat eonstpala oaah ofehop boeausa 
thsr aoofe Is a junction; tbs other Is about 
all fi-tuplos of segaent-lsbal . pairs whoro 
thj stroaa seguents aeot la a Junction and 
tho labels are possible for that type of 
junotloa. Let a stroaa junction be a ver- 
tex at whlob 3 stroaa segments aeet and at 
least ono segaent Is identified as visible 
rivers. Let J be tLe set of stroaa junc- 
tions, S ■ (t, ..., Ml be the set of stroaa 
segaents at stroaa junctions, . L • 
(upstream, downstreaal, and X be the set of 
juretlon patterns la Table 1. Let f: J - > 
X be the function that assigns labels to 
functions, and h: X -> L* be the funotion 
that specifies 3-tuples of labels that can 


meet at each type of Junotloa. For each a 
in X, we dtfme Tx and Is as follows. 

Ts e {(a. ,s.,s.) Is, ,a,,r t aeet la junc- 
tion type tl 1 * 3 * * 3 

IS e ( ( S j r )^ , Sj » lg • , 1 j ) | (if ,Sj,lj) 

0 fx and (l^lj.lj) « h (s)> 

Let T ■ U iM T f(g) and I e 0 -M »f(,)* 

T consists of all 3-tuplos of stroaa seg- 
aents that constrain each other, aad I is 
tho corresponding constraint relation. The 
labollag problea can be described by a ooa- 
patlblllty nodal (S, J , T, I), which is a 
particular lastaaee of the general consis- 
tent labollag problea tlarallck aad Shapi- 
ro, I $7 9 aad 1 9 AO ] • Because we believe 

that there are aany spatial inference prob- 
lea wbleh are instances of consistent la- 
beling problea, in the next section we de- 
scribe the fora of the general consistent 
labeling problea as given by ^llmana, Her- 
slick aad Shapiro £19*31. 


CiailBtiLt Lebellns mad Soafciol tenon . 
lna HbAAJL 

Let Q be a set of objects sailed 
units, aad L be a set of possible labels 
for those units. Let T a (flf a 0) be 
the collection of those subsets of units 
froa 0 that autually constrain oee anoth- 
er. That is, if f « (u. , u_, . . . , u k ) 
is an eleaeat of T, then not all possible 

labelings of u. u- are legal label- 

lags. Thus there is at least one label ea- 
st gone nt l , i_, , . . , L so that u. 
having labil if, u, having libel l,, . 

, u k having label l k is a forbidden 
labeling. T is called the unit oonstralat 
set. finally, let I £ Ig I giO t l, g 
single-valued, aad Ooa(g) 9 T) be the set 
of unit-label aapplngs in which constrained 
subsets of units are napped to their allo- 
wable subsets of labels. If g • 
( ( u I , 1. ) , (u, , 1-) , * * • , (u._,l k )l Is an 
element of s, then u f , u. , ... ; u. are 
distinct units, {u.,u,, f . . , u^} is an 
eleaeat of T aeaulng u., u. , . . . , u. 
autually constrain on* 'another, and uf 
having label 1., u_ having label i,, . . 

, and u k havlng^label l k are aj ' sim- 
ultaneously allowed. 

In the consistent labeling problea, we 
are looking for functloss that r.sslgn a la- 
bel in L to each unit i» ? and satisfy 
the constraints iaposed by T and I. That 
is, a consistent labeling is on* which when 
restricted to any unit constraint subset in 
T yields s sapping in I. In order to 
state this nor* preolsely, w* first define 
the restriction of a napping. Let h:0— >L 
be a function that ears eaoh unit in a to 
a label in L. Let f a 0 be a subset of 
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the units. The rsstrietlon h ! . (rend h 
restricted by f ) Is dsflnsd by hi* * 
{(u,l) « b ! u C fj. With this notation, 

vo dorms s consistent labeling os follows. 

1 function h:0 — > L Is s conslstont 
lsbollnf If snd only If for every f C T, 
h! f is on oloaont of 2. 

A« example la (Ivon bolow. Suppose tho 
Inputs to tbo probloa sro os follows: 

0 « Ct, 2. 3, 2, 5) 

L » (i, b, c) 

T « { {11, unory constraint 

(1, 21, binary constraints 

(2, 51, 

{1, 3, *11 tornary constraint 

fi « { {( 1 ,a] 1 , }( 1 ,M 1 , unary constraint 
C(1, a), (2. all, 

[(1 ,a) ,(2,b)l , 

{( I , b ) , ( 2, b ) } , binary constraints 

1(2, a), (5, all, 

((2, b), (5, oil, 

(( 1 ,*) ,{3 »a) , (b,c) 1 , ternary constraints 
{(1.b),(3.a),(t,a)/l 


tlao. 

Tbo spatial reasoning aodol (B , P, L, 
T, I, Ew) uan bo applied to fln'i tho flow 
directions of visible rivers. 0 oontslns 
tho units of visible rivers plus tho units 
of Invisible stroaas intersecting tho visi- 
ble rivers at Junotloca. P contains all 
tho properties dotootablo froa tho border 
segaents. L is (Opatroaa « 1 , Downstroaa a 
2}. T contains the junction relations. I 
contains the relations of legal flow direc- 
tions defined In Table 1. Sw is the nuabor 
of tlaes inconsistency occur* at Junctions 
normalized by the total nuaber of Junc- 
tions. 


£• Conclusion 

To detect stress network in LARDSAT, 
both visible rivers and invisible streaas 
need to be distinguished. Visible rivers 
can be detected by both spectral and 2D 
spatial lnforaatlon. Roweve’*, the detec- 
tion of Invisible streaas needs 3D spatial 
lnforaatlon. 


Then h a {(1,a) (2, a) (3, a) (t,e) (5, a)} 
is a consistent labeling. To see this note 
that h I jjj ■ (( t ,a) } , b ! ^ 2 } * 

{ { 1 , s) , t 2 , a) ) , h ! |2 jj » { (2, a) , (5 «a) } , 


S ii .Uuhv 


I. 


{( 1 ,a) , (3 »*) ,(*,e)l are 


If having (l t , .... 1^) applied to (u t , 
.... uj when (u., 1., .... u_, 1,) Is not 
in R is allowed with'a penalty, the process 
is called inexact consistent labeling 
[Shapiro and Harallck, 1981]. In order to 
Include these aapplngs, as error function 
Sw is defined. Let Ew : T X L fl -> [0, 1] 
be the error weighting function. Ew (u , 
..., u B , l , .... 1 } is the error which 
occurs when labels (l" .... l ) are ap- 
plied to (u^-r. . ••» Ugj . The napping h : 17 
-> L is an 1 inexact consistent labeling if 
for all (u. , ..., u.) in T, the suanations 

of Ew (u t , ..., Ug, h(u,j) h (u M )) is 

within sons upper bound. 6 

In spatial reasoning problems, aany 
spectral and geoaetrical properties can be 
detected for the locally detected units. 
Soae frequently used properties ire average 
gray level, size, and shape descriptors. 
Let P be the set of properties. The spa- 
tial reasoning aodel is (U, P, L, T, R, 
Sw). 0, L, T, Ew have the sane meanings as 
before; however, tne elements iu R now have 
-he fora ( u . , p . , 1 . , • • • , '!«, p *. , lu) whe- 
re p, is the list of property value ranges 
for All the properties in P for unit u., l 
* 1 to H. It Beans that if the property 
values of u, are within the ranges speci- 
fied by p^ for i « I to I, and (u,, ..., 
u.jj is contained in T, then it is legal to 
label lj to Uj , . .., lg to at the sane 


For invisible streaas, ridge end valley 
segaents aust first be detected and then an 
elevation growing aodel can be used to as- 
sign relative elevations to thaia. Interpo- 
lation oan generate surface elevation at 
all locations from the known values at 
ridge and valley segaents. Froa this ele- 
vation surface a valley network can be gen- 
erated easily. Another wny to fora a val- 
ley network is to create gap units as in 
Figure 7 during elevstlon growing. Local 
lnforaatlon including rules in Table 1 and 
other knowledge can be used to determine 
the Invisible streaa network as a subset of 
the valley network. The flow directions of 
invisible streaas cose directly froa the 
relative elevations of valley segaents. 

For the visible rivers part, the con- 
sistent labeling based spatial reasoning 
aodel can be used to find the flow direc- 
tions of visible rivers whose units are as- 
sumed to have all constant elevations in 
the elevation growing sodel. 

Based on the consistent labeling model, 
two types of spatial reasoning aodels can 
be formulated. If one is only Interested 
in the classification or labeling or the 
two dlaenslonal space so that roads, build- 
ings or other ground objects can be identi- 
fied, the aodel is called a 2D spatial rea- 
soning aodel. Cne exaaple is the aodel 
discussed in Section 8. If, in addition to 
the classification, soae feature values 
such as terrain elevations are needed over 
the two dlaenslonal spaoe, it is called a 
3D spatial reasoning aodel. In other 
words, in the 2D spatial reasoning aodel, 
the output specifications are symbolic; in 



tha ^0 spatial raaaoaiac aodal, the output 
specifications art nuaorlcal. 
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An application of the 20 spatial rea- 
soning nodal wbloh va hopa to report on 
soon la to reaognlse tha ground objects In 
an urban area by segaeatlag tha aarlal ia- 
aga Into regions, •aaaurlng tha properties 
of these reglons t fornulatlng constraints 
la T ana I, and applying the nodal. An ap- 
plication of the 30 spatial reasoning nodal 
uhloh ulll be dlsousaed la Haag's fortbooa- 
lag dissertation will be to find the beat 
set of segaent slopes so that the estlaated 
relative valley ridge elevations are as ac- 
curate as possible. 
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